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Abstract 

We re-examine a two-dimensional forest-fire model via Monte-Carlo simulations and show 
the existence of two length scales with different critical exponents associated with clusters 
and with the usual two-point correlation function of trees. We check resp. improve previ- 
ously obtained values for other critical exponents and perform a first investigation of the 
critical behaviour of the slowest relaxational mode. We also investigate the possibility of 
describing the critical point in terms of a distribution of the global density. We find that 
some qualitative features such as a temporal oscillation and a power law of the cluster-size 
distribution can nicely be obtained from such a model that discards the spatial structure. 
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1. Introduction 



Dynamical systems that are naturally close to or on their critical point have been proposed 
as an explanation of the appearance of power laws in nature [1,2]. This 'self-organized 
criticality' still consists to a large extent in the study of toy models and has many open 
or controversial questions. The lattice models can be grouped at least into three classes. 
A first class contains models with a local conservation law like the famous sandpile-model 
[1,2]. Recent experiments on piles of rice [3] partially confirm the general theoretical 
predictions by exhibiting power laws, but also show that the existence of power laws in real 
systems depends on microscopic details such as the aspect ratio of grains of rice. Models 
of evolution [4, 5] are one example in a second class where the dynamics is specified in 
terms of a globally selected extremal site. This extremal dynamics can be used to obtain 
some general statements for the complete class of models [6]. A model for self-organized 
criticality that is known as the 'forest-fire model' is a member of a third class of models 
that have parameters which can be tuned close to the critical point in a natural manner, 
but unlike in the two previous classes cannot be entirely discarded (In dealing with this 
model one should be aware that it is highly idealized and not expected to describe real 
forest fires). The precise version that we study in this paper has first been introduced 
in a short note [7] and arises as a certain limit of the more general model proposed later 
independently in [8]. 

The two-dimensional forest-fire model has been already discussed very controversially, 
mainly on the basis of Monte-Carlo simulations where usually the accuracy of the pre- 
dictions was the issue. An originally proposed version [9] did not show the desired critical 
behaviour [10, 11], and it was necessary to introduce lightnings [8]. Subsequently, Monte- 
Carlo simulations have several times lead to values for the critical exponents which had to 
be corrected later on [8, 12, 13, 14]. Here we add to this discussion by re-examining some 
of the quantities investigated only in one of the earlier works [12]. We were motivated by 
a study of the one-dimensional case [15] where we had discovered the existence of two dif- 
ferent length scales. The analogous question in two dimensions has been addressed in [12], 
but there it seemed that the two scales are proportional to each other. Here, we present 
more accurate simulations that demonstrate these two length scales to be different also in 
two dimensions. As by-products we also check or improve estimates for other exponents 
(see in particular [14]). A final subject is the approach to equilibrium which seems not 
to have been addressed in a similar way before. In a second part, we try to discard the 
spatial structure and introduce a global model similar to the one of [15] for one dimension. 
On the one hand, some of the qualitative features of the stationary state at the critical 
point (e.g. the power law of the cluster-size distribution) can nicely be described by such 
a global model. On the other hand, there are discrepancies in quantitative details and 
the range of such a simplified model is very limited in comparison to the full model. We 
believe that this is an important point, e.g. because it has been suggested in [16] that 
power laws in nature might arise from a global ('coherent') driving that does not see any 
spatial structure. Our findings here and in the one-dimensional case [15] demonstrate that 
the full model is not only different from but also richer than the simplified model. This 
is analogous to the result of [17] that versions of certain models of evolution with and 
without spatial structure lead to different results (at least if examined closely). 

We now define the model before we proceed with a presentation of our simulation results 
in the next Section. The forest-fire model is defined on a cubic lattice in d dimensions. 
Any site can have two states: It can either be empty or it can be occupied by a tree. The 
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dynamics of the model is specified by the following update rules (following [7, 12, 13, 14] 
In each Monte-Carlo step first choose an arbitrary site of the lattice. 

a) If it is empty, grow a tree there with probability p. 

b) If it is occupied by a tree, delete the entire geometric cluster of trees connected to it 
with probability /. This corresponds to a lightning stroke with subsequent spreading 
of the fire. 

A rescaling of the probabilities p and / just amounts to a rescaling of the time scale, and 
in particular leaves the stationary state invariant. We exploit this to set p = 1. There is 
a critical point at / jp = 0, but the parameter f jp is relevant and it is not legitimate to 
consider the forest- fire model precisely at this critical point (compare [18]). 

2. Simulation results in two dimensions 

In the following we consider the two-dimensional version of the forest-fire model on a 
quadratic lattice with periodic boundary conditions. The linear size of the lattice will be 
denoted by L. So the volume V is given by V = L 2 . We use a 'global' time scale, a unit 
of which is defined by the number of Monte-Carlo steps needed in order to visit each site 
on average once, i.e. a unit of global time consists of L 2 Monte-Carlo steps. 

In order to do the simulation efficiently also for large systems and small f /p one has to be 
careful and use e.g. bitmapping technologies. For details on the implementation compare 
the WWW page [19]. Using this program, the simulations of this paper took about four 
months of CPU time on 150MHz DEC alpha workstations. 

We have investigated mainly systems of linear size L = 16384 and parameter values 10~ 2 > 
f/p > 10 -4 . For a simulation, one random initial condition with density p = 1/2 was 
chosen. In order to equilibrate the system, it was left to evolving freely for at least 15 global 
time units. This equilibration time was increased to 25 global time units for / ' jp < 3- 10 -4 
and to 35 global time units for f/p = 1 • 10 -4 (these times were adjusted according to 
the observed time evolution of p, see also Section 2.3 below). After this, the system was 
iterated further for another 60 to 90 global time units (120 units for f/p = 3 • 10~ 4 ). 
During this period, measurements were made as global averages at intervals of usually one 
global time unit (for the density 200 times more frequently). This amounts to at least 
60L 2 measurements for each quantity of interest. Note that we use only a single run. 

2.1 Correlation functions 

Let us first discuss the simulation results for the usual tree correlation function {T(x)T(x + 
y)) (T(x) = 1 in a configuration with a tree at site x, T(x) = if the site x is empty). We 
have only investigated displacements y along the vertical axis (which can be treated partic- 
ularly efficiently using bitmaps). The treatment of the data resulting from the simulation 
is straightforward because it can nicely be fitted by 

C(y) := (T(x)T(x + ye 2 )) - (T(x)) 2 = ae'*'* (2.1) 

in a suitable interval y m i n < y < y max - The parameters a and £ were estimated by 
taking the logarithm of the r.h.s. of (2.1) and then performing a linear regression for 

x ) To be precise, the simulations in [13, 14] have been performed according to slightly 
different rules, because they aimed at investigating only quantities associated with clusters. 
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2/mm < V < !/max- These bounds are chosen such that the approximation (2.1) by a single 
exponential function is good. The lower cutoff can be chosen small (y m i n ~ 20) independent 
of f/p. An upper cutoff y max has to be imposed at distances of 4 to 7 times the correlation 
length £ because then statistical errors become large. 
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Fig. 1: The two correlation lengths £ and £ c . For £ c all estimates are 
on lattices with L = 16384 and the different symbols correspond to 
different ways of treating the data. 



Fig. 1 shows the results for £ obtained in this manner on lattices with L = 8192 (diamonds), 
L = 16384 (crosses) and L = 17408 (triangles). One observes that they are in good 
agreement with the form 

«~(£f T . (2.2) 

Performing linear regression fits on a doubly logarithmic scale one finds 2 ) 

v T = 0.541 ±0.004. (2.3) 

This is close to the result vt = 0.56 found in [12]. Comparing the values for £ obtained from 
simulations on different lattice sizes one can see that they may have a residual statistical 
error of around 2% which is about the same as the scattering of the individual data points 
around the line (2.2). Since also further data for larger f/p is consistent with the form 
(2.2) and the value (2.3), this result for ut and its error bound can be considered reliable. 

The estimates for the normalization constant a in (2.1) are compatible with a value a = 
0.030±0.001 independent of f/p. So, in the limit f/p — > the two-point function seems to 



l ) Error estimates are always the la confidence interval of a fit unless discussed explicitly. 
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tend to {T(x)) 2 +a for y > 100. In terms of the alternative ansatz C(y) = ay *? occ e y ^ used 

in [12] this corresponds to ?7o C c = 0, while the result found there was ?7 occ = 0.120 ± 0.015. 
We have also looked into the possibility of a power-law correction factor using our data 
for L = 16384 and f/p < 3 • 10 -4 . With the estimates for £ shown in Fig. 1 one finds 
that e y ^C(y) ~ y- - 11 for y < 20, i.e. for small y there is indeed a power-law correction 
factor with an exponent that is consistent with [12]. On the other hand, for 50 < y < 4£, 
the function e y ^C(y) is flat - its smallest values are around 0.027 and its maximal values 
around 0.031. This clearly contradicts a power-law correction factor with ?7 occ 0.11 in 
the large-distance asymptotics. Thus, for the large-distance behaviour ?7 occ = seems to 
be correct, while the power-law correction factor observed in [12] applies to small distances. 

We now look at a second quantity, namely the 'connected correlation function' (T(x)T(x + 
y)) c describing the probability to find two trees at positions x and x + y inside the same 
cluster. This quantity is usually referred to as the two-point function in the context of 
self-organized criticality and often is also the only correlation function that is investigated. 
We demonstrated in one spatial dimension [15] that the length scales associated to this 
correlation function and C(y) have different critical exponents while [12] suggested that in 
two dimensions length scales associated to different quantities are equivalent. This latter 
suggestion was based on simulations with lattice sizes up to 512 x 512 and f /p>5 ■ 10~ 4 . 

One of the main aims of the simulations presented here is to check if different length 
scales can be exhibited also in two dimensions after sufficiently improving the accuracy. 
As before, we restrict to displacements y along the vertical axis, i.e. we consider 

K(y):=(T(x)T(x + ye 2 )) c . (2.4) 

This quantity can be determined in the same run as the two-point function (T(x)T(x + y)), 
but it involves the additional effort of determining all clusters present in the system. 

Following [14] we associate a correlation length £ c to the second moment of K(y) via 

2 ._ Eg* V 2 K{y) 

The squares in Fig. 1 show values of £ c extracted from simulations with L = 16384 using 
this definition. One sees that these values are consistent with 

IV". («) 



and a (preliminary) value of 



0.5738 ±0.0013. (2.7) 



This agrees roughly with the value v = 0.60 found in [12], and also with the more accurate 
result in [14] which, including error bounds, is given by v = 0.580 ± 0.003 [20]. It should 
be noted that the lower bound y = 1 in (2.5) is crucial. Starting the summation instead 
e.g. at y = 10, one finds £ c m 24 at f/p = 10~ 2 and £ c m 219 for f/p = 10~ 4 instead of 
£ c 14 and £ c ~ 203, respectively. This in turn would make the value for v smaller. We 
will argue soon that the start of summation y = 1 is indeed the correct choice, but the 
impact of a modification here could give rise to a slightly larger error than the la interval 
for the fit given in (2.7). 
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We now proceed with a more detailed discussion of the form of K(y) which will also justify 
the definition (2.5). Eq. (2.5) is based on the following expected form of K(y): 

K(y) = a c y-^ e~ v ^ . (2.8) 

In particular, if e y ^ c K{y) agrees well with a power law, the determination of £ c via (2.5) 
is justified. Using the values of £ c given by the boxes in Fig. 1 one finds that e y /^ K{y) 
does indeed agree well with a c y~ v for y between 1 and several times £ c 3 ). The crosses 
in Fig. 2 show these estimates for rj. They obviously depend on f/p contrary to what one 
would naively expect. This points to systematic errors in the determination of rj which 
one may expect to become less important for smaller f/p where the power law is better 
visible. This suggests to make a 'scaling ansatz' rj(f/p) = rj + a (f /p) b to determine the 
value of 7] in the limit f/p — > 0. A least-squares fit gives the dotted line in Fig. 2. One 
finds a critical rj = 0.374 ± 0.011. 
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Fig. 2: Values for the exponent r\ obtained in two different ways 
(crosses and diamonds) and scaling fits (lines). 



Alternatively, one can fit all three parameters in (2.8) directly from the data. The values 
for i] obtained with this approach are shown in Fig. 2 by diamonds. A scaling analysis 
(the line with long dashes in Fig. 2) yields a critical ?] = 0.342 ± 0.008. Comparing this 
with the value found earlier, one finds a small disagreement. So, the error bounds of these 
two estimates for the critical r\ may be a little too optimistic because they do not include 
systematic errors. In order to be careful we give a final result that includes both estimates 
for rj and the two error ranges: 

3 ) The values of a c are all compatible with the //^-independent value a c = 0.20 ± 0.01. 
Therefore, at the critical point the power law K(y) = a c y~ v is expected to be valid for all 

y- 
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71 = 0.36 ±0.03. 



(2.9) 



This is compatible with the result rj = 0.411 ± 0.02 of [12] (although the error estimate 
of [12] seems a little too optimistic, possibly because it does not take systematic errors 
into account). One also obtains different estimates for £ c from the same fits that gave 
the alternative values for r\. These values are shown in Fig. 1 by 'x'. One sees that 
they are again compatible with the form (2.6), and obtains an alternative estimate for the 
associated critical exponent: v = 0.5777 ±0.0013. There is a minor difference between this 
value and the earlier estimate (2.7) indicating that we have indeed neglected systematic 
errors. Therefore we give a final result 



which includes the two direct estimates and their error bound. This final result is in 
excellent agreement with [14]. 

Comparing (2.10) with (2.3), we see that v ^ i>t, so that £ c and £ are basically different 
lengths. A difference v—vt ~ 0.04 was already observed in [12], but attributed to numerical 
errors because it was unexpected. 

It should be mentioned that our results (2.9) and (2.10) for r\ and v do not satisfy the scaling 
relation 2 — r\ = 1/v [12]. We can only speculate about the reason for the disagreement. 
For example, in the derivation of this scaling relation one assumes that r\ and v do not 
depend on the direction of the displacement vector y, and we have not checked whether 
this is indeed true. Another possibility is that we have still overlooked systematic errors. 
Inserting v (which is the more reliable value) according to (2.10) into the scaling relation 
2 — r\ = 1/v yields r\ 0.28 which is possible if our error estimate in (2.9) is by a factor of 
about 3 too small. 

2.2 Cluster-size distribution 

The distribution n(s) of clusters with size s arises as a by-product of the determination of 
K(y) during the simulations. We extract an exponent r from it following the lines described 
in detail in [12,13]. One introduces the quantity P(s) = J2 s '>s s'n(s'). Assuming that 
it behaves as P(s) = as 2_T e~ s / Smax , one can use three-parameter fits to obtain estimates 
for r and s max - Extrapolation of the values obtained in this manner from the simulations 
with L = 16384 yields r = 2.1595 ± 0.0045 for f/p -> 0. The values of r for f/p > 
approach this limiting value from above. An alternative way to extract a value of r is 
to look directly at n(s) and assume that n(s) ~ s~ T for intermediate s. Applying this 
second method to the same data and to some results for L = 8192 we obtain the estimate 
t = 2.159 ± 0.006 at the critical point. This limit is now approached from below and is 
in excellent agreement with the one obtained before. To be on the safe side we retain the 
value with the larger error bound as the final result: 



This value agrees within error bounds with most previous results [12, 21, 13, 14], but our 
error bound is considerably smaller than in most of them. 

When one extracts estimates for r from P(s) one also obtains estimates for s max . These 
values are in good agreement with the form s max ~ (/ /p)~ x with A = 1.171 ± 0.004. To 



v = 0.576 ±0.003 



(2.10) 



t = 2.159 ±0.006. 



(2.11) 
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check the validity of our error estimate we can use the scaling relation A = 1/(3 — r) [14]. 
After inserting (2.11) one finds the prediction A = 1.189 ± 0.008. This prediction agrees 
roughly with the direct estimate. However, there is a small discrepancy indicating that 
the estimate 

A = 1.17 ±0.02 (2.12) 

is more realistic. Within error bounds we find good agreement with the values of [14] and 
[12]. From the results (2.10) and (2.12) we obtain the fractal dimension p = 2.03 ± 0.04 
using the scaling relation p = \jv [14]. This does not quite agree with [14] where p = 1.96± 
0.01 was found, but would favour instead the expectation p = 2 of [13]. Unfortunately, 
with our simulations we had not aimed at determining p, and we are therefore not able to 
clarify this interesting point. 

2.3 Density and time evolution 

Here we study the critical behaviour of the stationary density as well as the temporal 
behaviour of the density. The latter yields the lowest gap in the spectrum of the time- 
evolution operator which is given by a straightforward generalization of eq. (3.6) in [15]. 
Some high relaxational modes of this time-evolution operator can be written down ex- 
plicitly as we demonstrate in Appendix A. In one dimension, the low-lying spectrum of 
this operator could be studied numerically [15], but this is not feasible in higher dimen- 
sions. Fortunately, simulations of pit) exhibit much clearer features in two dimensions 
than was the case in one dimension. This is illustrated by Fig. 3 which shows the initial 
time evolution of pit) after the system was started at p(0) = 1/2. 
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Fig. 3: The density p(t) as a function of global time t. Shown are 
results from simulations with f/p = 10 -2 (bottom, full line), f/p = 
10~ 3 (line with long dashes, middle) and f /p = 10~ 4 (top line with 
short dashes). 
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First we examine the critical behaviour of the stationary density of trees p(oo). The values 
in Table 1 are obtained by averaging p(t) over times t after the equilibration time, i.e. at 
or beyond the right border of Fig. 3. 
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Table 1: Estimates with L = 16384 for the stationary density of trees p(oo) and 
the lowest decay mode in p(t) of the two-dimensional forest-fire model. 

These values are in good agreement with 

Pc-p(oo)~^0 , 

where the critical density p c and the critical exponent 5 are given by 

p c = 0.40844 ±0.00011, (2.14) 
1/5 = 0.466 ± 0.004 . (2.15) 

These results agree within error bounds with the values obtained in [14] , but our bounds 
are smaller. In particular we can now rule out that 1/5 = 1/2 as proposed by [13]. 

Having determined p(oo), it is straightforward to extract the oscillation period and decay 
time of the slowest decay mode from p(t). First, one determines those t where pit) crosses 
the value p(oo). From the distances between these crossings the oscillation period T osc can 
be determined easily. Averaging 10 to 15 half oscillation periods estimated in this manner 
for a suitable interval of time in Fig. 3 leads to the values given in Table 1. One observes 
that the oscillation period equals T osc = 0.88 within error bounds for all f/p, i.e. the 
slowest relaxational mode oscillates with a constant frequency. This is to be contrasted 
with the one-dimensional case where simulations of p(t) clearly demonstrated that the 
oscillation period depends on f /p [15]. 

Finally, we extract the leading decay time T from p(t) according to the following proce- 
dure. At times t precisely in the middle between two subsequent crossings used for the 
determination of the oscillation period, the value of \p(t) — p(oo) \ is determined. One finds 
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(2.13) 



for these (approximately ten) values of t that \p(t) — p{po)\ ~ exp(— t/T) from which it is 
straightforward to obtain the estimates for T presented in Table 1 4 ). These values for T 
are compatible with a critical behaviour 

T~(£) (2.16) 

where the critical exponent ( is determined to be 

C = 0.314 ±0.013. (2.17) 

This exponent probably is another new exponent that is not related to the ones determined 
so far [12,14]. 

2.4 Summary of simulations 

Table 2 summarizes our results for the critical exponents of the two-dimensional forest- 
fire model. It also includes the critical density p c and the global oscillation period which 
seems to be independent of f/p. For comparison we have also included the results for 
one dimension [15,22]. In that case, the oscillation period diverges and we have listed 
the exponent rather than a period in Table 2. Similarly, the amplitude a vanishes in one 
dimension for f/p — > but is roughly constant in two dimensions. 
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exponent ~ 0.194 
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0.030 ±0.001 
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Pc 


0.40844 ±0.00011 


1 



Table 2: Summary of our results for the critical behaviour of the 
forest-fire model. 



4 ) These estimates also verify that our equilibration times are long enough, since we 
have equilibrated the system for at least 6T before starting to collect data. So, the non- 
stationary modes are damped by factors of at least exp(— 6) ~ 2 ■ 10 -3 during data collec- 
tion. 
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3. Global models 



We now wish to investigate through simplified variants of the model to what extent one 
can describe the stationary properties of the two-dimensional forest-fire model by global 
variables in a way similar to [15]. There it was shown that in order to describe the one- 
dimensional critical stationary state it suffices to know the relative weight of the sum of all 
configurations with a fixed number of occupied (or empty) sites. Thus one is dealing with 
a kind of grand canonical ideal lattice gas. Such a model has no intrinsic spatial structure 
and leads to a two-point function independent of the distance. It can therefore be used 
to describe the asymptotic behaviour of the critical C{y) of Section 2. Working with the 
(global) density of trees p, one has to specify the probability distribution p(p) and obtains 
C(y) = (p 2 ) ~ (p) 2 f° r V m the thermodynamic limit. This is positive for all continuous 
distributions. 

In two dimensions the limits f/p — > and L — > oo do not commute with each other (in 
contrast to one dimension). Therefore, perturbation theory cannot be used to compute 
p(p), and in fact we do not know of a good analytic method to determine p(p) from the 
rules described in the Introduction. Therefore we use heuristic arguments and simulations 
to discuss it. 

The forest-fire model reminds one of site percolation. Thus, it is natural to try to relate 
the forest-fire model to percolation and gain some insight from that. Attempts in this 
direction have been made e.g. in [8, 21]. We will also try to use some results of percolation 
theory, but we will follow a different route. Namely, we try to interpret the stationary state 
of the forest-fire model at the critical point as a suitable ensemble of percolation problems 
with a distribution p(p). 

It is known from percolation theory [23] that in a homogeneous configuration with a density 
p above the percolation threshold p perc two arbitrary sites are connected with a finite 
probability. In such configurations, the dynamics of the forest-fire model with arbitrarily 
small f/p > would very quickly destroy all percolating clusters and thus drive the density 
below the percolation threshold 5 ). Thus, the probability p(p) to have a global density 
above the percolation threshold p perc must vanish in the forest-fire model, i.e. p(p) = for 
P > Pperc- In one dimension one has p per c = 1 and in two dimensions p per c = 0.592746 [23]. 
In all simulations p(t) was well below this percolation threshold at all times (compare Fig. 
3). 

It is useful to visualize the stationary state of the forest-fire model in order to gain some 
intuition. Fig. 4 shows an area of 760 x 472 sites at t = 66 during a simulation with 
L = 16384 and f/p = 10" 4 (compare also Fig. 2 of [24], Fig. 1 of [14] and Fig. 6 of [13]). 
One observes that at a certain fixed time the system consists of rather well-defined patches 
with different mean density of trees. These patches are not to be confused with single 
forest clusters, they usually contain many such clusters. Their typical size increases as f/p 
becomes smaller, which reflects the divergence of correlation lengths. Looking at the time 
evolution of such a state 6 ), one observes an increase of density due to growth of trees that 
is constant throughout the patches and that lightning strikes essentially only the patches 

5 ) This implies that the density p cannot be continuous at f/p = for d > 1 and is the 
reason why perturbation theory cannot be used in higher dimensions. 

6 ) On XI 1 platforms, such a visualization is possible with the code used for the simu- 
lations presented in this paper. This code is available on the WWW [19]. 
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with the highest density. After a lightning has struck such a patch, a new patch with a low 
density is created. This new density is not really zero because there are always some trees 
in the patch that are not connected to the cluster which is destroyed. The idea now is that 
the important information about the critical stationary state is given by the distribution 
p(p) of densities in these patches and that nothing essential changes if we replace such 
patchy systems with an ensemble of systems of global density p occurring with the same 
probability p(p). Of course, it remains to be tested to what extent this picture works. 




Fig. 4: Snapshot of an area with 760 x 472 sites in the stationary state 
during a simulation with L = 16384 and f/p= 10 -4 . Trees are black 
and empty places white. 



3.1 A simple model 

Let us make a very simple model based on these ideas. Assume that below the percolation 
threshold p per c trees just grow (with probability p = 1) and no lightning strikes. As soon as 
the percolation threshold is exceeded, lightning strikes immediately and destroys all trees 
in the system (not just those in the percolating cluster). In order to determine p(p) we first 
compute the mean lifetime of a configuration with density p, assuming that no lightning 
strikes (this consideration will also be useful later on). Such a configuration lives precisely 
n local updates if n — 1 times an occupied place and then an empty place are selected. 
Because of the above assumptions there is no correlation and therefore the probability of 
this to happen is given by p n_1 (l — p). This yields the expectation value t(p) of the lifetime 
as 

oo 1 

t{p) = Y,np n -\l-p) = —. (3.1) 

n=l ' 

If no lightning strikes, the probability p(p) to find a configuration with density p is pro- 
portional to the lifetime of a state with this density. Taking into account that lightning 
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strikes all trees at p perc , we find for this simple model 

^ u j P ^ Pperc 7 

with the normalization constant given by M~ x = J Pporc dp (1 — In one dimension, 

(3.2) agrees with eq. (5.4) of [15] which was derived there using a different argument. 

This simple model has the following periodic time evolution: Trees grow until the density 
reaches the percolation threshold. Then lightning strikes and the process is restarted with 
a completely empty system. The average time needed for such a cycle is precisely the 
global oscillation time and is given by J Pperc dp t(p). Inserting t(p) according to (3.1) and 
the value of p peT c in two dimensions yields an oscillation period of T osc ~ 0.898 which is 
in very good agreement with what we found in simulations (compare Table 1). The mean 
density is given by f Q dp pp(p) from which one finds a critical density p c = 0.340 . . . (this 
deviates notably from the result (2.14) found by simulations of the full model). Finally, 
the probability to find two trees at arbitrary places is given by the second moment of 
p(p), i.e. by f Q dp p 2 p(p). So, the two-point function exceeds the value p\ by an amount 
a = 0.0289 . . . which agrees within error bounds with what we found by simulations for 
the large-distance asymptotics of the two-point function. Cluster-type quantities are not 
accessible as easily and would e.g. require again Monte-Carlo simulations. 

Although this simple model yields very good values for two quantities, its failure to give 
the correct p c is not surprising. Firstly, we have neglected the fact that lightning can 
already strike configurations with p < p per c even for arbitrarily small f/p > (compare 
Appendix A). Secondly, lightning does not really lead to the completely empty system, 
but usually leaves some isolated trees or small clusters in the patch behind. Unfortunately 
we do not know how to treat either effect analytically. This lack of knowledge also has the 
effect that we can in general not compute an oscillation time from p(p) although we do of 
course still think of the system as evolving in cycles (compare also section 6.3.2 of [25]). 

3.2 Realistic distributions 

Next we try to obtain a realistic p(p) from Monte-Carlo simulations. Measurements of the 
global density cannot be used to extract p(p) from a simulation, because p fluctuates only 
very little around its mean value for sufficiently large systems (see Fig. 3). Therefore, one 
has to look at the distribution of local densities. In a large system, areas with different 
local densities coexist and we are interested precisely in these local fluctuations and not 
just the global average. We have decided to divide the system into 16 x 16 plaquettes and 
use the distribution of the average density per plaquette. This size of the plaquettes was 
chosen because then their linear extent is much smaller than the correlation lengths and 
they contain sufficiently many sites to obtain a fairly smooth distribution. Fig. 5 shows 
a result p re (p) obtained in this manner using the parameter values closest to the critical 
point, namely / jp = 10~ 4 and L = 16384. Samples were taken at the same 90 times where 
also the correlation functions were determined, amounting to a total of almost 10 8 samples 
for local densities. The normalization in Fig. 5 is such that Ylr^oPre( r '/256) = 1. 

As explained above, the first and second moments of p(p) are related to the mean density 
and the asymptotic value of the two-point function. From this one finds 

p(oo) = 0.4044 ... , a = 0.031.... (3.3) 
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The way we have determined p Te (p) ensures that the first moment indeed equals the value 
obtained by directly taking a global average. The slight difference between (3.3) and the 
value for p(oo) in Table 1 is due to the fact that the latter is based on a much larger 
amount of configurations. As for the simple model presented before, the prediction for 
C(y) = a ~ 0.031 agrees within error bounds with what we expect for the asymptotics of 
the two-point function at the critical point. 




0.2 0.4 0.6 0.8 1 

P 

Fig. 5: Result for the distribution p re (p) on 16 x 16 plaquettes in a 
simulation with L = 16384 and / jp = 10 -4 . 



Even though this basic test yields good values, one should be aware that examining the 
system only in windows blurs the distribution in Fig. 5 in several ways. Firstly, looking 
through a window containing only 256 sites, one finds a smearing of the density by Ap pa 
0.02 just because of statistical effects. Secondly, such windows may accidentally intersect 
the boundary between two patches with low and high densities. This has the effect that 
p(p) is estimated too large for intermediate p and too small for the extremal ones. Our 
choice of 16 x 16 plaquettes is designed to make both effects reasonably small and is about 
the best we can do. 

Let us now discuss Fig. 5 keeping these effects in mind. For not too small values of p, 
Pre(p) increases slowly and reaches a maximum around p ~ 0.54. Around p ~ 0.62 there 
is a sharp decrease. The broad distribution shows that fluctuations of p are important. A 
peak just below the percolation threshold p pe rc and a steep decrease above it correspond 
to our expectation. However, there is still a substantial contribution to p re (p) above p pe rc 
which is not explained by windowing effects. This is due to the patchy structure of the 
system: Finite patches with p > p peTC are not destroyed instantly, but rather live for a 
time which is the longer the smaller these patches actually are. One could also say that the 
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patchy structure demonstrates that the system is actually correlated. In particular at small 
distances (j/<20), the two-point function C{y) retains a y-depence (compare Section 2.1) 

for f/p — > and thereby exceeds the asymptotic constant a. This correlation is necessary 
to observe a non-trivial distribution of local densities, but also leads to contributions to 
p(p) above the percolation threshold. 

So, if we want to work with a 'realistic' p(p), the best we can do is to proceed with the one 
shown in Fig. 5. Alternatively, one can work with an approximation to this distribution 
where one suppresses the undesired p(p) above the percolation threshold by hand. One 
such approximation we have examined is a linear one, i.e pn n (p) ~ p for p < 0.59 and 
Piin(p) = for p > 0.59. This linear distribution yields p c 0.395 and a ps 0.019 (both 
somewhat too small). 



1 — ,— n — 1 — i— n — 1 — i— n — 1 — i— n — 1 — i— n — r 




s 



Fig. 6: The cluster-size distribution obtained from a simulation of the 
full model with L = 16384 and f/p= 10~ 4 (full line), the one based on 
Pre(p) (long dashes) and the result obtained frompii n (p) (short dashes). 

The determination of the cluster distribution n(s) and K(y) is a complicated combinatorial 
problem which we solve again using simulations. One creates configurations with densities 
distributed according to the given p(p) and then measures n(s) and K(y). We have created 
100000 configurations on a 512 x 512 lattice distributed according to p re (p) and 70000 
configurations on a 1024 x 1024 lattice for pii n (p)- Fig. 6 shows the cluster-size distribution 
n(s) obtained in this manner together with the one obtained from a simulation of the full 
model. For small cluster sizes (s < 100), all three distributions are close to each other. 
However, at larger s the distributions based on a globally given p(p) decay faster than the 
true n(s). The corresponding exponent is r ~ 2.48 for the distribution p Te (p) and r ~ 2.44 
for pim(p) - both much closer to the mean-field value r = 5/2 [21] than to the true value 
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(2.11). In Fig. 6 one also observes a peak in the cluster-size distribution corresponding to 
p r e(p) for 1 • 10 5 < s < 2.5 ■ 10 5 ~ 512 2 , i.e. just below the volume of the system. This is 
due to the non- vanishing of p re (p) for p > p pevc which leads to clusters spanning a finite 
(and large) fraction of the system. 




Fig. 7: The full line shows the correlation function K(y) obtained 
from a simulation of the full model with L = 16384 and / jp = 10 -4 . 
The result for p Te {p) is shown by the line with long dashes, the one 
obtained from pr in (p) is indicated by the shorter dashes. 



Fig. 7 shows the probability K(y) to find two trees at distance y inside the same cluster. 
Here, the results obtained from the two f»(p)'s deviate more notably from the result ob- 
tained by simulation of the full model (full line). As noted above, the distribution p Te (p) 
gives rise to percolating clusters and produces a constant background in K(y) of approxi- 
mately 0.0552. After subtracting this background, one can fit K(y) with (2.8) for £ c ps 180 
and r\ ps 0.92. This value for rj is more than twice as large the true one (2.9). For pn n (p) 
one finds good agreement with the form (2.8) for £ c = L/2 and rj ~ 0.95. 

3.3 Summary and outlook on global models 

We have shown that the distributions p re and p\\ n lead to a power law for n(s). One can 
check that the same is true for other p(p)'s that are cut off at p peTC in a way similar to 
piin- Thus, a power law in n(s) arises automatically from a description in terms of global 
quantities and need not be a signal for criticality in a conventional sense. However, in all 
examples for p(p) discussed so far, we obtained values for r and r\ that are unsatisfactorily 
larger than those of the full model. Nevertheless, a description in terms of p(p) can still 
be forced to work because for two-dimensional critical percolation one has r ps 2.055 and 
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rj 0.208 [23] - values which are smaller than the ones in Table 2. So, one can obtain 
the desired value e.g. of r by peaking the distribution p(p) more prominently just below 
Pperc- Adjusting just r to its correct value can be expected to also give a reasonably 
good approximation for rj. Afterwards, both p c and a may be tuned to the desired values 
by adding another peak in p(p) for smaller densities (which does not contribute to large 
clusters and thus affect the asymptotics of cluster quantities). However, there does not 
seem to be a natural way to make these adjustments. In particular, local densities above 
Pperc do exist in the full model which would have to be discarded by hand in a global model 
in order to obtain the correct K(y). 

In one dimension the spatial structure becomes irrelevant at the critical point [15]. We 
have seen in this Section that this can be generalized to the qualitative features of the 
critical correlations in two dimensions, but not to the quantitative details. In contrast to 
the one-dimensional case we had no analytical tools at our disposal and have therefore not 
been able to derive the distribution p(p) of local densities explicitly. The introduction of 
block-spin variables is reminiscent of real-space renormalization group ideas and it would be 
interesting to see if they can be used to find p(p) . However, one would have to go beyond the 
block-spin renormalization-group study of [18]. Firstly, one would have to admit densities 
different from or 1 for the block-spin variables, and moreover the dynamics should not 
be treated just in mean-field approximation. 

Two-dimensional percolation is believed to be conformally invariant (see e.g. [26]). The 
globalized models are just ensembles of percolation problems and should therefore be con- 
formally invariant as well. It would be interesting to know if also the stationary state 
of the full model in two space dimensions is conformally invariant, even if the standard 
techniques of conformal field theory would probably not say much about quantities like 
cluster sizes. 

4. Conclusions 

In this paper we have again looked at several aspects of the two-dimensional forest-fire 
model. Firstly, we have shown that the two length scales £ and £ c have different criti- 
cal exponents. That this might be possible had been suggested by a study of the one- 
dimensional model [15] which illustrates that one-dimensional systems can provide useful 
insights because of their relative simplicity even if one is actually interested in higher- 
dimensional versions. This result shows that in generic non-equilibrium systems geometric 
objects and the usual (occupancy) correlation functions can behave completely differently. 
For an equilibrium system as the Ising model such an observation was already made some 
time ago in [27]. In this case, percolation occurs away from the critical point, i.e. the two 
length scales are so different that they diverge at different temperatures (see e.g. [28, 23]). 

In order to show that vt ^ v we had to improve the error bounds of earlier investigations. 
As a by-product we have also improved the accuracy of other critical exponents. It may be 
possible that one could still improve the error bounds by another digit using optimized code 
on today's most powerful computers. Historically, Monte-Carlo simulations have already 
several times lead to values for the critical exponents that had to be corrected later on 
[8, 12, 13, 14] and as we have shown here, some of them were still not treated adequately. 
Therefore, it may be desirable to perform yet another independent verification of the results 
presented here, but a further increase of accuracy may not be necessary for this end. 

The second part of the paper focussed on a globalized model. We found that one can 
easily obtain power laws in the cluster-size distribution by discarding the spatial structure, 
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i.e. by making the usual two-point function independent of the spatial coordinates. This 
generalizes a result obtained for one dimension in a previous paper [15] and is line with the 
observation in [16], based on a different one-dimensional model, that one can obtain power 
laws in clusters or avalanches by global ('coherent') driving. However, some quantitative 
predictions of the globalized model did not work out satisfactorily. One reason is that the 
full two-dimensional forest-fire model has non-trivial two-point functions at least at small 
distances which is also reflected by the existence of patches with local densities above the 
percolation threshold. In addition, the full model exhibits many critical exponents that 
cannot be described by a global model. 

A study of the usual correlation functions would also be desirable in other models of self- 
organized criticality where they have not yet been investigated. This could help to clarify to 
what extent the process of self-organization can be regarded as a global phenomenon. Two- 
point correlation functions would also be important quantities to examine in experiments. 
It would e.g. be interesting to extract the spatial correlation functions of the local heights 
and slopes from the experimental data of [3] . 

Even for the two-dimensional forest-fire model there are still many issues we have not 
looked at, including e.g. finite-size effects. We have only looked at the regime fV ^> p 
which is close to the thermodynamic limit. One could also look at a different limit, namely 
fV <^ p where the first-order approximation of [15] applies independent of the spatial 
dimension. In this limit, trees grow until the lattice is full and after a certain time of rest, 
lightning destroys all these trees and the process starts again. It would be interesting to 
see how this behaviour crosses over to the critical behaviour studied here as the volume of 
the system is increased, and if finite-size scaling can be observed. 
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Appendix A. Some exact relaxational modes 

The time-evolution operator for the forest-fire model in general dimensions d can be written 
down along the lines of Section 3 of [15]. Then one can obtain exact excited states using 
the same argument as at the end of Section 3 of [15] (cf. in particular eq. (3.15) loc. cit.). 

The crucial step is to consider a class of configurations that contain a single cluster of 
trees and where each empty place is a neighbour of a tree. The last condition means that 
the configuration still consists of a single cluster if a tree is grown at any of the empty 
places. Now we consider states that are built out of such a configuration but have non- 
zero momentum P^O. Lightning strokes in such configurations lead to the completely 

— * — * 

empty system. In states with P^O the completely empty system is reached by lightning 
with coefficients that contain a sum over all roots of unity such that lightning maps states 
constructed in this manner to zero. Thus, the only terms in the image come from growing 
a single tree in a configuration. Now, as in eq. (3.15) of [15] one can write down eigenvalue 
equations that are upper triangular in the number of empty places N. The diagonal terms 
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in these equations are already the eigenvalues 



A N = Np + f/p(V-N). (A.l) 

One must have N > 1 because the completely full system cannot be used to build a state 
with non-zero momentum. 

The number of empty places N is restricted by the above conditions. For d = 1 only the 
two configurations considered in [15] meet the requirement of consisting of a single cluster 
and remaining in this class after growing an arbitrary tree. Thus, N < 2 for d = 1, and 
the density p of these two configurations tends to 1 in the thermodynamic limit. 

For general d note first that the connectedness of the cluster implies that the trees must 
form at least one-dimensional objects. This is most efficiently implemented by grouping 
the trees in straight lines. In the hyperplane perpendicular to these lines, a tree can have 
2d — 2 empty places as neighbours. Thus, at least one place out of 2d — 1 must be occupied 
by trees in order to meet the requirements. This yields N <V(1 — l/(2d— 1)). Equivalently, 

all configurations containing only a single cluster that retain this property after growing a 
tree must have p > l/(2d — 1). For d<3we can indeed specify configurations that saturate 
this lower bound. In d = 1 this limit is p = 1 and corresponds to the thermodynamic limit 
of an almost full system. In d = 2 a configuration with p = 1/3 is given by lines of trees 
that are mutually separated by two lines of empty places. The lower bound p = 1/5 in 
d = 3 is saturated by arranging the lines in such a way that their mutual position in the 
perpendicular plane corresponds to the moves of a pawn in the game of chess. Of course, 
the lines have to be connected by a d— 1 dimensional object, but this does not change the 
value of p in the thermodynamic limit. 
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